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Abstract 

A high-ranking goal of interdisciplinary modeling approaches in sci- 
ence and engineering are quantitative prediction of system dynamics and 
model based optimization. Quantitative modeling has to be closely re- 
lated to experimental investigations if the model is supposed to be used 
for mechanistic analysis and model predictions. Typically, before an ap- 
propriate model of an experimental system is found different hypothetical 
models might be reasonable and consistent with previous knowledge and 
available data. The parameters of the models up to an estimated con- 
fidence region are generally not known a priori. Therefore one has to 
incorporate possible parameter configurations of different models into a 
model discrimination algorithm which leads to the need for robustifica- 
tion. In this article we present a numerical algorithm which calculates a 
design of experiments allowing optimal discrimination of different hypo- 
thetic candidate models of a given dynamical system for the most inap- 
propriate (worst case) parameter configurations within a parameter range. 
The design comprises initial values, system perturbations and the optimal 
placement of measurement time points, the number of measurements as 
well as the time points are subject to design. The statistical discrimina- 
tion criterion is worked out rigorously for these settings, a derivation from 
the Kullback-Leibler divergence as optimization objective is presented for 
the case of discontinuous Heaviside-functions modeling the measurement 
decision which are replaced by continuous approximation during the op- 
timization procedure. The resulting problem can be classified as a semi- 
infinite optimization problem which we solve in an outer approximations 
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approach stabilized by a suggested homotopy strategy whose efficiency is 
demonstrated. We present the theoretical framework, algorithmic realiza- 
tion and numerical results. 

1 Introduction 

High-ranking goals of interdisciplinary modeling approaches in the natural sci- 
ences are quantitative prediction of system dynamics and model based optimiza- 
tion. In particular in modern systems biology a related issue is to link molecular 
attributes to dynamic mechanisms and functional properties at the system level 
in order to mechanistically understand emerging functionality. For these pur- 
poses, mathematical modeling, numerical simulation and scientific computing 
techniques are indispensable. Quantitative modeling closely combined with ex- 
perimental investigations is required if the model is supposed to be used for 
sound mechanistic analysis and model predictions. 

Typically, before an appropriate model of a system is found different hy- 
pothetical models might be reasonable and consistent with previous knowledge 
and available data. The goal of this article is to derive, develop, implement 
and apply a numerical algorithm which calculates in a suitable sense an optimal 
design of experiments which allows the best discrimination of different hypo- 
thetic candidate models in form of ordinary differential equations (ODE). The 
algorithmic idea is to iteratively separate the response of different models by 
use of variations of experimental conditions and perturbations to the system. 

To discriminate a set of candidate models against a given set of experimental 
data likelihood ratio tests based on bootstrap methods have been described in 
the literature, see e.g. [52], [15] or [55]. Ranking methods like Stewart's method 
([M]) or the well known Akaike information criterion (see e.g. [IS]) are popular 
as well in the field of biological modeling. Applications can be found for example 
in gl] or HD]. 

In contrast to these approaches our work deals with the problem of designing 
experiments so that statistical methods can be exploited in an optimal sense for 
model discrimination. This differs from the approach to find an experimental 
design to best estimate the parameters of a model for a given experimental 
system in terms of criteria characterizing the confidence regions [231 [7] E] . 

Different approaches to design experiments for model discrimination exist. 
Besides optimization methods (see e.g [29], [19], [47] or [26]) a model-based 
feedback controller see e.g [3] and Markov chain Monte Carlo sampling methods 
|34j have been used to construct an appropriate design. An overview of various 
experimental design techniques can be found in [2"7] . 

Here, we present a robust numerical optimization algorithm which calculates the 
optimal design of experiments allowing the best discrimination of different candi- 
date ODE models. An appropriate model and its parameters up to an estimated 
confidence region are not known a priori. Therefore one has to incorporate pos- 
sible parameter configurations of different models into a model discrimination 
algorithm. The aim is to calculate the most discriminable response of different 
models for the most inappropriate parameter configurations within a parameter 
range via a worst case estimate. In that context inappropriate parameter con- 
figurations refer to the case when different models have calibrated parameter 
values such that the models have the most similar response. This worst case 
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estimate leads to the formulation of a maxmin optimization problem. Building 
on our previous work |43j we present an algorithm to compute robust optimal 
experimental designs. For the robustification we set up an outer approximation 
approach stabilized by a homotopy strategy. 

The article is organized as follows. In Section 12.11 we give a brief overview of 
so called Kullback-Leibler(KL)-optimality as discussed by Lopez-Fidalgo ct al. 
[52"] in the context of model discrimination. In Section |2~21 we derive our optimal 
design criterion by use of KL-optimality. In Section [3] the theoretical framework 
for the calculation of a robust design via solution of a maxmin optimization 
problem is presented. The numerical implementation is discussed in Section 
13.11 The homotopy solution strategy is presented in Section 13.21 We demon- 
strate applications of the algorithmic framework for two test cases from biology, 
an allosteric metabolic enzyme model for glycolytic oscillations and a model 
describing signal sensing in dictyostelium discoideum, results are presented in 
Section HI 

2 Statistical Basis of Model Discrimination 
2.1 KL-optimal design 

In this section a model discrimination criterion based on the Kullback-Leibler 
(KL) divergence called KL-optimality as discussed by Lopez-Fidalgo et al. [3J] 
is introduced. Lopez-Fidalgo et al. [32] demonstrate that KL-optimality is con- 
sistent with T-optimality [5] and generalized T-optimality [35] which are well 
known model discrimination criteria based on statistical tests. 

We introduce the concept of a probability space and formally define the KL- 
divergence. 

Definition 1. The probability space is a triple (fl, T, P) consisting of 

• a non-empty set (sample space), 

• a a -algebra T C 'P(fi), E <E T is called an event 

• a probability measure P : T — > [0,1]. 

Definition 2. Two probability spaces J 7 , Pi), i — 1,2, are called absolutely 
continuous with respect to each other, in symbols Pi = P 2 , if J E £ T : 
(Pi(E) = AND P 2 (E) ^ 0) OR (Pi(E) ^ AND P 2 {E) = 0). 

The Radon-Nikodym Theorem allows a representation of a probability mea- 
sure via a measurable probability density function. 

Theorem 1. (Radon-Nikodym) 

Let X be a probability measure such that A = P\ , A = P 2 - Then X-measurable 
functions fi : — > M, i — 1, 2, called generalized probability densities, exist 
which are unique up to sets of measure zero and non-negative, such that 

Pi(E)= f fi(x)d\(x), 1 = 1,2, (1) 

J E 

for all E G T . 
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A proof of this theorem can be found e.g in [T2"| . 



In the following we use X for the generic variable and x for a specific value 
of X. If Hi, i = 1, 2 is the hypothesis that X is from the statistical population 
with probability measure Pi, the mean information for discrimination in favor 
of Hi against H 2 given x £ E € J ' , for P\ is given by the Kullback-Leibler 
divergence. 

Definition 3. The Kullback-Leibler (KL) divergence is given by 

0, ifP 1 (E) = 0, 

with 

dPr(x) =h(x)dX(x). (3) 

When E is the entire sample space £1, we shorten the notation to I(P± : P2). 
For discrete sets E the integral is substituted by a sum. 

For details we refer to [28] . 

Now assume that the sample space is split into two disjoint sets E\ and 
E2, = £1 U £2. We define a statistical test procedure to choose between 
hypotheses H\ and Hi by accepting H\ if x € E\ and accepting H2 if x G E^. 
Assuming that one of the hypotheses has to be true we treat H2 as the null 
hypothesis and call E\ the critical region. The following wrong test decisions 
can occur. 

Definition 4. Incorrectly accepting H\ although H2 is true is called the type I 
error. The probability that this error occurs is given by 

a = Prob(:r £ E l \ H 2 ) = P 2 {E l ) . (4) 

Definition 5. Incorrectly accepting H2 although H\ is true is called the type II 
error. The probability that this error occurs is given by 

p = PTob(xeE 2 \H 1 )=P 1 (E 2 ). (5) 

Wc assume that the test is repeated n-timcs and denote by O n a sample of 
n independent observations. 0\ represents a sample of a single observation. j3 n 
is defined as the corresponding probability of an error of type II which depends 
on the number of independent observations and the splitting of the probability 
space O into disjoint sets E\ and E2- 

The following theorem demonstrates an asymptotic relation between the KL- 
divergence and the minimum possible probability /3* of an error of type II with 
respect to all possible splittings E1UE2 = ft with given a = Prob(x € -Ei | i?2 ) = 

p 2 (e 1 ) ng. 



4 



Theorem 2. For any value of a, say <xq, < ao < 1, 

lim (ftf'" = e-x&.P.M (6) 

n— >oo 

A proof of this theorem is given in [T51 [2H] ■ 

Assuming probability models for the outcome of a data measurement ex- 
periment depending on experimental design parameters £ G S C M d , this the- 
orem justifies the KL-divergence to be an appropriate objective functional for 
model-based computation of an optimal experimental design for discrimination 
between model hypotheses. For a design with the largest possible value of X the 
asymptotic probability /3* of encountering an error of type II becomes minimal 
with respect to all possible splittings E\ U E<i = with given oto- We indicate 
the dependency of the KL-divcrgcncc on the design by Z(P 2 : P\ , 0\ ; £) . Our 
aim is to derive an algorithm to calculate the optimal design £ G S such that 

1 = arg maxZ(P 2 : P 1: Or, £)• (7) 

An extension of the case to test a simple null hypothesis against a simple alter- 
native hypothesis to the more general case of both hypotheses being composite is 
generally of interest. This includes the situation to test whether given measure- 
ment data can be explained best by one out of a finite set of probability models 
based on measures P ri , r\ G {1, ...,Mi} parametrized by parameters 9 ri G ri 
where n C M Pri is the set of all possible parameter values to parametrize 
P ri and Mi is the cardinality of the set of probability models, against the hy- 
pothesis that the measurement can best be explained by another one out of a 
second finite set of probability models based on measures P r2 , parametrized by 
parameters 8 r2 G 6,- 2 where 6 r2 C K p '' 2 and r 2 G {1, M 2 }. 
By calculating 

£ = ar g max min min T(P r2 (9 r2 ) : P r , (6 ri ), @\\ (8) 
?es ne{i,...,M 1 } e^ee^ 
r 2 e{i,...,M 2 } e r2 ee r2 

we can get a robust worst case estimate of an optimally discriminating design 
for the case of composite null and alternative hypothesis. 

In practical applications a simple strategy to sort different probability models 
into null and alternative hypothesis would be to first rank all models according 
to the existing measurements and then set the best model as null hypothesis 
and the others as alternative hypothesis. The development of a suitable and 
efficient strategy is subject to further work. 

It should be noted that the presented criteria is not symmetric, i.e. the null 
hypothesis is favored. In the case that both hypotheses might be equally rea- 
sonable we suggest to use the symmetrized version of the KL-divergence, i.e. 

as optimization objective instead. 
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2.2 Derivation of the optimal experimental design crite- 
rion 

In this section we derive a numerically computable optimization objective func- 
tional based on the framework of KL divergence. The derivation is motivated 
by the requirements of biological in vitro time series experiments modeled by 
kinetic ODE systems. In most situations such experiments are time and cost 
consuming. Therefore a central issue is to get the most information out of a 
single time series data measurement experiment taking place within a given 
fixed time span [0,T cn d]. This means that in an optimal experimental design 
the most informative measurement time points for one measurement run have 
to be calculated in such a way that only one measurement at one time point can 
be performed. Often, an experiment cannot produce measurements in a time 
continuous way. Therefore we assume that there has to be a minimal time span 
AT for the separation of subsequent measurement time points. This contrasts 
to the usual approach to associate weights to a discrete or continuous time de- 
sign scheme, see e.g. [5]. Additionally, the initial species concentrations of the 
participating species should be chosen in a most discriminating way. 
A commonly used experimental practice is to combine kinetic time series mea- 
surements with perturbation stimuli like external adding of species quantities. 
From the model discrimination point of view the optimal time point of pertur- 
bations and the optimal species quantities to be added should be determined. 
We further assume that a measurement cannot be done at the same time as a 
perturbation. 

In the following we translate these experimental conditions into a statistical 
model. Given the measurement time-vector t £ K™ with entries U for the n 
measurement time points U,i £ {l,...,n} such that ti + ± > ti, the "internal" 
model response vectors y\. at measurement time U for the parametrized prob- 
ability measures P rj , rj £ {!,..., Mj}, j £ {1,2} are given by 



Vlj ■= Vr, (U-i,k, 1 + Ci-i,O rj ), (10) 
1 + c. 



where y Tj (ti-i , ti , y l r . 1 + Ci_i , 8 rj ) £ R mr j are the solutions of the initial value 
problems 

/^(y^AJ, t€[ti,4 (ii) 



dt 

with initial state y rj (t l ~ 1 ) := J/£ + c;_i at end time ti where to := and 
Co := 0. The vectors a £ l™ m " denote species quantities the system can be 
perturbed with at time points ti where m max is the maximum dimension of the 
models, i.e. 

m max := max max m r .. (12) 

je{l,2}re{l,...,M,} 1 

f* hs (-, ■) are the right hand side functions of the ODE models, yj £ I™ m " 
denotes the initial species concentration vector of the entire experiment which 
is for all models the same, i.e. y®. :— yi. By m m i n we denote the minimal 
dimension of the models, i.e 

m ra in := min min m r - (13) 

iS{l,2}re{l Mj} 3 

We do not assume that m n 

the "redundant" entries in Cj and yj are "ignored" 
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To each ODE model we associate an observable function O rj : M" lr j — > K° 
which describes an experimental observation explained by that model where o 
denotes the dimension of the experimental observation. The expected observa- 
tion O l r . of the rj-th model at time point i is given by 

(>: ■■ o, ;//; >■ (w) 

Let 0' i denote an observation at measurement time point ti. By assuming 
that the measurements at successive time points U are independent with nor- 
mally distributed error vectors e* . G R" 1 with zero mean and variance functions 
v rj (0\. , ti, 9 rj ) 2 we get for the regression models 

O 4 * = Tj + e r . (15) 
the model probability densities f rj (-', ■) at measurement time point £j given by 

/-,«V,<*,) = i—e-Wr^rw-o*^ (16) 



with |w lirj | := rife=i v r- (O r - > **> ), where u^. (O* . , £*, ) denotes the fe-th en- 
try of the square root of the variance functions v r . (O*. , U, 9 r . ) 2 , and diagonal 
matrices V? G M oxo with diagonal entries [V^] := (O* , ij, ri )) 2 



We generally allow different error models. The error models might dependent 
on the expected observations O r . , the time U and possibly on parameters Q Tj . 

For the sake of notational simplicity we define 

/. ■<>' ■ ■ f, (0 U ;O r )■ (17) 

For a full measurement run containing n measurement time points we get the 
probability density models 

n 

f, (O): Hi'. >(>' ;■ (18) 

However, by assuming such a model probability distribution we still allow that 
two measurements arc separated by a time span less than AT. 

To overcome this problem we extend the probability spaces f2j = K° of a mea- 
surement at one measurement time point by one-element-containing sets A/", 
to 

fi< = MM (19) 

where £l{ is the disjoint union of Oj and Mi- The element of the set Ni with 
measure P(M%) G [0, 1] represents the event "no measurement", i.e. O ti G Afi <=>■ 
"no measurement performed at time point ti" . 

In order to derive measures on fij, i = 1 , . . . , n that allow for a density func- 
tion representation according to the Radon-Nikodym theorem (Theorem [I}, we 
introduce the Heaviside-function 

H : R+ — > [0, 1] (20) 
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with 

By use of this Hcavisidc-function and er-algebras T%, where T contains the 
Lebesgue measurable sets on fli and additionally the union of these with the 
set Afi, we define probability spaces (fli,Ti, Pi. rj ) with measures 

Pi, Tj -Ei&Ti^ Pi, rj (Ei) e [0, 1] . (22) 

Three cases have to be distinguished: 1. Ei C £1;, 2. Ei C Afi, 3. Ei n Oj 7^ 

and EinAfi^Q. 

For case one with c f2i we set 

hr ] (E l ):='H(U) f f r .{O u )dO u . (23) 

For case two with Ei C A/^ we set 

Pi, ri (Ei):=l-U(ti)- (24) 
For case three with Ei n 0, 7^ and S< n A/i 7^ we set 

A,r,(#i) := H(ti) f f T] {O u )dO u + (1 - H(ti)) • (25) 

By introducing these modifications the probability models based on measures 
Pi <r . do not depend on measurements which are performed in less than AT time 
after the previous measurement any more. 

To take into account that a species concentration perturbation can only be 
applied if no measurement is done at the same time, the same procedure is 
repeated with the additional Heaviside-function 

«(«>-{; 

The measures Pi <rj are defined in the same way as above replacing %(ti) by 
U[ti)M{ci). 

For specific r\ and inserting the probability models Pi >ri and respectively 
Pi,r 2 m to the KL-divergence (Definition [3]) using Xi := Pi_ r2 an d the additivity 
of the KL divergence for independent events one gets the following expression 



l(P r2 :P ri ,0 1 )=J2 



i=i 



H(*i)W(ci)/„(0**) log j M^iS^^l 



(i-nmici)). ^ 1 -^^ 



(27) 

where c„ := 0. With log(l) = this simplifies to 



l(P r2 :P ri ,0 1 ) = Y,U{t i )n{a) J /r2 (0^).log||^j|dO*'. 



(28) 
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By inserting the normal distribution (JTBJ) in (|28[) one gets 
n 

l{P r2 :P n ,0i)=^K(ti)W( Ci )- 



j=i 



/r 2 log 



i e - h (oi 2 ~o^) T v; 2 (oi 2 -oH ) 



27r|i) iir2 | 



J e -i(0^-0'l) T V^(0* 1 -0*i) 



This is equivalent to 



i(p r2 :P ri ,o 1 ) = J2^(u)n( Cl ) 



v k ri (Ol. 1 ,U,9 ri ) 

■['~ V* 2 (0}. 2 ,t h 9 r2 ) 



Al 



k > 



dO u . 

(29) 

(30) 



with 



Al :=- 



1 

'(^(o^tiA,)) 2 
l 



«(0* t.AJ) 2 vl ^J* 



L^j:-2^ ri j fe L^j fe +L^j: 



(31) 



and where [0\ k gives the k — th entry of the observation vector O. A\ reduces 
using the well known moments of the normal distribution to 



Ai = i 



( \o l r2 \ I - 2 [oi 2 \ I + [oi 2 \ I + « (oi 2 , u, e 2 )y 

( l°k\ I - 2 [o i ri \ k L°rJ fc + L°rJ I + « (o; 2 , «*, ^ )) 2 ) 

(32) 



This further simplifies to 



(lpij 2 fc - 2 [oi i \ k [oi 2 \ k + [oy 2 fc + «(o; 2 ,t t A 2 )) 2 ) 
1+ (<(o* ^AJ) 2 



(33) 



Substituting A\ back into (|3Q|) we get 

n o / 

l(P r2 : P ri ,Oi) = Y, ^feM*) E lo S 



i=l 



fc=l 



/ <(o; 1 ,t i ,g r2 ) N 

«* (O* ,*i,0 ri ) 



-i 



(34) 
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This reduces to 



1 ™ 

l(P r2 :Pn,0i)=-5^(ii)tt(ci). 



i=l 



E 



«(o^ 2 ,^,^ 2 )) 2 + (Lo; 2 J fc -KJj' 



-2 log 



(35) 



This criterion has to be maximized with respect to the initial concentration 
vector yi, the measurement time point vector t and the system perturbation 
vector c, thus £ := (yi,t,c) € 5 C M d . 

For our optimal experimental design we generally start with a large number of 
measurement time points. By use of the Heaviside functions the number of mea- 
surement time points is reduced such that for ti — < AT the corresponding 
measurement time point is "turned off" . 

These Heaviside-functions %{•) and %{■) can be replaced by any appropriate 
continuously differentiable switching functions with range space [0, 1]. 
It should be noted that we assume that we have the same time discretization 
for measurements of different species and the addition of further species quan- 
tities. This assumption is practical especially for the application to in vitro 
experiments performed by biologists. For introducing arbitrary generic con- 
trols we need a more general formulation of time schemes, i.e. simultaneously 
time schemes which are independent of each other. One is associated with the 
controls, others may be associated with distinct observables which might be 
measured indepently. The incorporation to the presented framework is subject 
of further work. 



3 Solution of the max-min optimization prob- 
lem 



We formally state now the experimental design optimization problem Pe: 

max t (36) 

(T,e)escR i + 1 

subject to 



min I(P r2 A,) :P ri (0 ri ),O i; O-r>O, n e {1, Mi}, r 2 G {1, M 2 }, 
e ri ee ri 

n 

At, = Tend, 

yf n <yi< y? iax , 
o < At < < max , 
o < c < c max , 
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withe := (yi,M,c) G K d , At t := ti-Vi and 9 := {(6 ri , e r2 )|r 3 - e {1, M 3 }}. 
The auxiliary variable r £ 1 is used to transform the maxmin optimization 
problem ([8jl to a maximization problem with an infinite number of inequality 
constraints. The remaining constraints model the feasible range of experimental 
setups. 

To solve optimization problem (|36[) numerically by applying efficient derivative 
based algorithms we replace the Heaviside functions H(U) and H(ci) in ([55)1 
by continuously differentiable approximations, parametrized hyperbolic tangent 
functions of the form 



tanh( 



6(Atj-bi) ■ 



1 ~ tanh(- 
— and % (ci) = 



6(cj-fc 2 ) ' 



1 



(37) 



The parameters ai.2 characterize the width of the transition region between 
and 1. The parameters 61.2 determine the center of the transition region (see 
Figure [T]) . By setting the parameters in an adequate way arbitrarily close ap- 
proximations of the Heaviside functions can be generated. A different approach 
to handle the discontinuous Heaviside functions would be to introduce binary 
variables and treat the resulting problem as Mixed Integer Nonlinear Program- 
ming problem. The drawback of this approach is that its solution can become 
very expensive. There seems to be little theoretical work in literature on Mixed 
Integer maxmin problems and an efficient solution strategy is not obvious in 
that case. 




width a - 




Figure 1: Switching functions: the left switching function is used to guarantee 
that only one measurement is done at a time point, the right one is used to 
guarantee that if a perturbation is done at a time point no measurement is done 
at the same time point. 



In literature, problems as (|36[) fall into the class of semi-infinite inequality and 
equality constrained optimization problems (SIECP) [551 . 
Several methods to solve such SIECP problems are available, an overview can 
be found in J5TJ [55] . We choose the method of outer approximation [?TJ [55J , 
whose origin can be traced back to cutting plane methods for convex problems 
|38j . This approach is beneficial in the presence of a complex inner problem, in 
our case the robustification against the parameters (9 ri and 9 r2 . The outer ap- 
proximation algorithm iteratively solves discretized finite counterparts P§n of 
the semi-infinite problem Pe in each step N, successively refining the discretiza- 
tion Q N until a sufficient approximation of the original problem Pe is reached. 
For problem (|36[) this means that in each iteration of the outer approximation 
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scheme ri and 0,. 2 are replaced by finite approximations 8^ and respectively 
6^ with N := {(6,^,0^)1^ G {1, ...,Mj}}. This relation between the semi- 
infinite problem and an infinite sequence of finite problems can be formalized 
in the theory of consistent approximations and epi-convcrgence 35, 36, 37 , 381. 

We use a modified version of Algorithm 3.6.4 in [38] where "Step 1.", the cal- 
culation of augmenting vectors 0^ +1 and 0^ +1 to construct 

:={ClU^ and 9f 2 +1 := U 6* (38) 

is realized by 

: = arg min X(P r2 (6r 2 ) ■ PnVnlOnbr), (39) 

with ^tv denoting a locally optimal design of the previous problem Pg« ■ The 
algorithmic scheme is as follows: 

Algorithm 1. 



• Data. Choose (o£S and a sequence {cat}^ =1 with e^r > and e^r 4- 0- 

• Step 0. Set N = 1, set 9°. := 0. 

• Step 1. Calculate <d N according to \38\) and \39\) . 

• Step 2. Calculate approximate solution of Pgw such that 

-e N , (40) 

and the equality and inequality constraints in problem Pgw are fulfilled up 
to ejy. 

• Step 3. Replace N by N + I, and goto Step 1. 



iPg (•) : W l+1 — > R<o denotes the optimality function associated to problem 
Pgw, see Theorem 2.2.24 in [38]. The optimality function ^g (•) is always non 
positive and directly related to the first order generalized Karush-Kuhn- Tucker 
(KKT) conditions, i.e. 'J'g ((r, £)) = if evaluated at a generalized KKT point, 
see Theorem 2.2.19 in [ggj^ 

Assuming that l(P r2 (9 r2 ) : P ri (0 ri ), X ; £) and V(I{P r2 {e r2 ) : P ri (0 ri ),0i;£) 
are Lipschitz continuous on bounded sets with respect to £ and 8 rj and <d rj arc 
compact any accumulation point of Algorithm [1] fulfills the generalized KKT 
conditions, compare Theorem 3.6.5 in |38j . 

To calculate 6* +1 in Step 1. of Algorithm Q] we use on heuristic base a simple 
random search approach coupled to a local optimization method, i.e. we have 
randomly generated P different start values in Q rj , from which we have started 
the local optimization method for parameter estimation. The best value out 
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of the P trials is chosen to augment the set 9^. Of course there are more 
sophisticated approaches to search for a global minimum for a review see e.g. 
[U, but at this point an effective calculation of Step 1. of Algorithm [1] was 
not our primary goal. For the local parameter optimization we use the same 
optimization method as for Step 2. in Algorithm [TJ 

In our implementation we use a fixed e at the desired final accuracy to solve 
problem Pg N in every loop of Step 2. of Algorithm [TJ i.e. ej^ — e, N > 0. 
In that way Step 1 . of Algorithm [TJ gives a worst case estimate of the KL di- 
vergence I(P r2 (0 r2 ) : P ri (0 ri ),0i;£jv) for the current design £jv with respect 
to 6 rj . Therefore for practical application the algorithm can be stopped if the 
worst case estimate of KL divergence for the current design £jv is big enough 
although no local optima might be achieved during optimization. As stopping 
criterion of Algorithm [TJ we use: 

Algorithmic Stop Criterion. 

Stop after Step 1. of Algorithm]]^ if 



S> min min l(P r2 {9 r2 ) : P ri (9 ri ), Or, 6v-i) 

ne{l,...,Mi}0 ge"- 1 

^•••' M2 V, e e£- 



(41) 

mm mm l(P n (6> r , ) : P ri (9 ri ), £>i; Cw-i) =: A_r G , 
nG{i,...,Mi} o ri ee ri 
r-y |i \h \ e r2 ee r2 

where 5 is a small positive constant, then consider (#^ +1 , 9^ 2 +1 ), T\ £ {1, Mi}, 
r% £ {1, A/2} and^N as (approximate) solutions of problem Pe, else goto Step 
2. and calculate a new design £,n+i- 

This stop criterion has also been used in [T^l [33] . We call the distance Ana 
given by (TTTj) . robustification gap. 



3.1 Numerical solution of the problem Pgiv 

We have implemented the resulting optimization problem in a multiple shooting 
setup (see for example [T5J[TTJ[T3]). In our multiple shooting setup the whole in- 
tegration interval [0, T en( j] is subdivided into several subintervals by introducing 
auxiliary multiple shooting node variables s rjt i_i, j £ {1,2}, rj £ {1, Mj}, i £ 
{1, n}, I £ {1, A}, on each of which an independent initial value problem 
is solved. Each end point of a subinterval corresponds to one measurement time 
point. Matching conditions which enter the optimization problem as additional 
equality constraints assure continuity of the state trajectory from one subinter- 
val to the next. 

To incorporate the perturbations c matching conditions 

Sr :i ,i,l-yr :i (ti-l,ti,Sr j ,i-l,l,d' r:l ) =0, (42) 

s rj ,i,i, denoting the multiple shooting nodes with s rj fl.i = Hi are modified to 

s Tj ,i,l -y r AU-i,U, s jt i-i t i,6 l ) = Ci, i £ {1, n — 1}, 

~l ( 43 ) 

Srj,n,l Urj (j*n— 1 j j Sj,n— 1,1 j ^r? ) ^* 
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A graphical scheme of the multiple shooting setup is shown in Figure [5] 
Instead of evaluating the objective functional (f3"5"j) by use of the values y\. , given 
by the solution of the initial value problem (|lip , depending on the parameters 
rj , (f3"5j) is evaluated by use of the auxiliary multiple shooting node variables 
s rj ,i,i replacing the values y\. with s Tj ^\ respectively. The dependency of ((35)) 
on s r .ij is indicated by I{s r - L ,-,i, s r2 ,-,i)- 



modified 
matching condition 



Fixed end time T 




time 



Figure 2: Scheme of the multiple shooting setup for computing the experimental 
design to discriminate two models. A dot denotes one measurement time point. 
The black solid line corresponds to model 1 and the gray dashed one to model 
2. 



The overall optimization problem can be stated as 



max t 

r,yj,A£,c,s 



subject to 



~~sr 



frf S {y,O l rj ), te[t<_l,ti], VriiM-l) :=*r„i-l,J 
&rj,i,l IJrj {pi— 1 j ^ij $rj ,4—1,2) ^Tj ) ' ) 
$rj ,n,l Vvj ipn— 1 j ^71 j ^rj ,n— 1,/j ) ^' 

Sr f ,o,/ = 2//, 
2/} nin < yi < yf ax , 
< At < i max , 
< c < c max , 
s min , < s r . ; i < s max ,, 

I(*n,.,J,*r a ,.,0-T>0. ri G{l,...,Afi}, r 2 s{l,..,M 2 }, 
with je {1,2}, r, G {1, Mj}, ie{l,..,n}, J e {1, ...,JV}. 



(44) 



(45) 



We have implemented this problem within the interior point optimization pack- 
age IPOPT (SOJElE], using the linear solver MA27 [Mj. Usually for the solution 
of the KKT system within the direct multiple shooting ansatz the linear system 
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is solved utilizing tailored structure-exploitation, e.g., condensing or high rank 
updates. See for example [3D]. Since speed aspects are not our primary concern 
we rely on the sparse solver MA27 instead of developing a tailored solver for 
this problem class at the current stage. 

All derivatives up to second order, which are used for the calculations of the 
Hessian needed for a robust performance of IPOPT, are calculated by automatic 
differentiation using CppAD, [pJ|S]. 

For the solution of the differential equations within the optimization problem, 
which are commonly stiff in chemical and biochemical applications, we have 
implemented a fully variable step, variable order (order 1 to 6), Backward dif- 
ferentiation formulae (BDF) method, based on Nordsiek array polynomial inter- 
polation similar to the EPISODE BDF method by Byrne and Hindmarsh |16j . 
but with the step size selection strategy of Calvo and Randez [TT| . 
For the generation of sensitivities we have adopted the sophisticated principles 
of internal numerical differentiation developed by Albersmeyer and Bock (5J [T] 
in forward and adjoint mode. 

The idea of this principle is instead of calculating the sensitivities by use of 
the sensitivity differential equation, to directly differentiate the BDF integra- 
tion scheme by automatic differentiation, which we implemented using CppAD 

mm- 

According to some notes in the PhD thesis of Albersmeyer [1] we also have im- 
plemented the possibility to control the step size scheme not only by the local 
truncation error of the nominal trajectory but as well by the local truncation 
error of the sensitivities generated by the forward mode of automatic differen- 
tiation with respect to the sensitivity differential equation, which has shown by 
numerical experience to improve the robustness of the optimization approach. 
A different approach would be to use collocation, i.e. to incorporate a full dis- 
cretization of the ODEs into the optimization problem, see e.g. [TT]. Since the 
kinetic ODE systems in the focus of our applications are usually stiff, we prefer 
adaptive time integration. 



3.2 Stabilizing homotopy method for subsequent Pgjv+i 

By solving the subsequent optimization problems Pq N +i with an interior point 
code like IPOPT [50j [51] initialized with primal and dual variables of the pre- 
vious problem or with primal variables only, one often observes that the new 
solution may differ significantly from the previous. This is due to the fact that 
the solution of the previous problem Pgw is infeasible for Pq N +i and thus the 
algorithm tries to find a feasible state before it proceeds to find a new optimum. 
This behavior is not desired in the context of an outer approximation algorithm, 
because convergence of the algorithm may be slowed down significantly. This 
circumstance originates from a jumping between vicinities of distinct local max- 
ima of problem (|36p . The discretization Q N of the robustification space may 
not be equally adequate for different local maxima. To overcome this problem 
we have implemented a heuristic homotopy method to gradually introduce the 
additional constraints 

flV-i,r 2 ( T j£)-N"+l : = Z(Sr u -,N+l, S r2 ,-,N+l) - T > 0, (46) 
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of problem Pgw+i- We replace ffn,r 2 ( T :£)jv+i by 

5n,ra(T,C; k) n+1 :=X(s ri) .,jv+i,s ri!l . 1 jv+i) - T + (1 -«)/>> 0, (47) 

with homotopy parameter kg [0,1] and p is a constant which has to be set such 
that g ri ,r 2 { T -i £; K )./v+i are inactive for re = at the initial design £jv- We choose 



K is a save guard factor we set empirically to K = 1.4, which worked well 
in practice for our examples. For re = the augmented optimization problem 
should be easily solvable within a few iterations by performing a warm start from 
the solution of the previous problem. By increasing the homotopy parameter 
to k = 1 the additional constraint is gradually introduced, which leads to a 
sequence of easily solvable subproblcms whose solutions stay in the vicinity of 
the solution of the previous problem Pq„ . A similar homotopy strategy can be 
found e.g. in |4T)] . 

4 Numerical results 

We have applied the algorithm developed in Section [5] and [3] to two example 
problems for which we present results in the following section. For the purpose 
of illustration we restrict ourself to the case that each hypothesis comprise only 
one model whereby we assume that only the alternative hypothesis is composite. 
We also assume that each species is "directly" measurable. We treat model 1 
as null hypothesis and model 2 as alternative hypothesis. 

4.1 Discriminating design for two models describing gly- 
colytic oscillations 

In the first test case for model discrimination we implemented the following 
models for glycolytic oscillations as described in [20] . 

Model 1 is an allosteric enzyme model with positive feedback under cooper- 
ativity and linear product sink. The differential equations for model 1 are given 



p to be 



p := K ■ max 

ri6{l,.,Mi) 
r 2 e{l,...,A/ 2 } 



min I(P ra (0 ra ) :P ri (e ri ),Oi;^) 




(48) 



by 




— = gi£70(a!i,7i) - fc s 7i, 
ai(l + ai)(l + 71) 2 



</>(<*!, 7i) 



il + (l + Ql) 2 (l+ 7 l) 2 ' 
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Model 2 is an allosteric model with positive feedback in the absence of co- 
operativity and the product sink is represented by Michaelis-Menten kinetics. 
The differential equations for Model 2 are given by 

_ = u - 0(0:2,72), 
— = g 2 0(o:2,72) 



(02,72 J 



dt [i + 72 ' 

0:2(1 +72) 



L 2 + (l + a 2 )(l+72)' 



ai ; 2 denotes the species concentration of the substrate 71,2 that of the product. 
For both models the inflow parameter v is the same and fixed to the value 
v = 0.22. It represents the inflow of substrate to the experimental system, a 
CSTR (continuously stirred tank reactor). 

The parameters a, q±, k s and L\ of model 1 are regarded as known. Their values 
are given in Table [TJ the parameters qi , t s , (1 and L2 of model 2 are regarded 
as unknown and subject to robustification. For the permitted parameter range 
see Table [TJ 



Model 1 


Model 2 


a 
0.92 


9i 

2.01 


k s 

0.11 


Li 
17206.10 


92 

[io- 7 ,ioo] 


[io- 7 ,ioo] 


I 1 

[10- 7 , 100] 


£2 

[100, 300] 



Table 1: Parameter values for the glycolytic oscillation models. 

For simplicity we consider the homoscedastic case with equal variances, i.e. 
vi = V2 = cr 2 . In this case I (Pi : P 2 ,Oi) reduces to, 

n 

X(Pi : P 2 ,Oi) =Y,U'(t i )U'(c i ) ((a\ - o^f + (7? - itf) . (49) 

i=l 

For this test case the homotopy strategy as presented in Section 13.21 is only ap- 
plied if the robustification gap Aug < 0.1, then the successive problem Pg N +i 
is calculated by use of the homotopy strategy with 30 homotopy steps, i.e. 
Kh = ft/30, h G {1, ...,30}. Otherwise the problem Pg N+ i is solved without ho- 
motopy strategy. For each subsequent problem Pq N +i the solution of problem 
Pqn is used as initial guess. 

We first present a robust design without the possibility to disturb the system 
by adding species at later time points. 

The design is calculated within a fixed time window i.e. T on <3 = 400. 100 equally 
spaced possible measurement points are defined in the initial state of the opti- 
mization procedure, the distance vector At between the time points is subject 
to design and each entry is restricted to Ati G [10 -7 , 10 19 ], i G {1, 100}. The 
disturbance vectors Cj are set to Cj = 0, i G {1, ...,99} and are fixed to model 
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the fact that no species disturbance is allowed. 

The initial species concentrations which are also subject to experimental design 
are restricted to aj £ [10~ 7 , 25] and 77 G [10 -7 , 25]. The initial values were set 
to ai = 15 and 7/ = 2. The parameters of the switching functions H'(ti) are 
chosen as a\ = 20.0 and bi = 10.0. The parameters of the switching functions 
U'(ci) are chosen as a 2 = 0.05 and 6 2 = 0.025. The algorithmic settings are 
summarized in Tabled 



Optimization settings 


Integrator settings 


P 
5 


6 

10~ 6 


IPOPT-tol: Step l./Step 2. 
10 _io /10 _8 


rclTol/absTol 

io- 12 /io- 12 


relTolScns/absTolSens 

io- 12 /io- 12 



Table 2: On the left hand side the optimization settings are listed comprising 
the IPOPT stopping tolerances for Step 1. and Step 2. of Algorithm Q] and on 
the right hand side the integration tolerances for the nominal trajectory and the 
first order sensitivities. We use the IPOPT option "honor_original_bounds=no" 
for Step 1. and Step 2. of Algorithm [1] 

A plot of the functions ai, «2 and 71, 72 in the initial state and for the solution 
of problem Pgj are shown in Figure [3] A plot for the same functions with the 
same solution design as for problem Pgj after the next robustification step is 
shown in Figure U] The final design is also shown in Figure 21 A plot of the 




50 100 150 200 250 300 350 400 50 100 150 200 250 300 350 400 

time time 



Figure 3: The model functions ax,Q!2 and 71,72 are shown before the opti- 
mization procedure (left) and after the optimization procedure of problem Pgj 
(right) for the glycolytic design setup without the possibility to disturb the 
system. One square represents one measurement time point. 

robustification gap A rq and as well for the objective value of problem Pgw for 
each iteration N of Algorithm [1] are shown in Figure [SJ A selection of design 
variables as solutions of problem Pgw is shown in Figure IHKlcft). 

In a second scenario we additionally allow for species perturbations. In this 
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Figure 4: The model functions ai, 0,1 and 71, 72 are shown for the same solution 
design as for problem Pg! after the next robustification step (left) and for the 
final design (right) for the glycolytic design setup without the possibility to 
disturb the system. One square represents one measurement time point. 



new scenario at the 21th 41th, 61th and 81th measurement time points, the 
system can get disturbed by additional species quantities. The free vectors Ci, 
i G {21,41,61,81} are constrained by Cj G [10~ 7 , 10]. The initial values are set 
to Cj = 1. The remaining conditions are as before, however we change the time 
vector bound constraints for i G {1,6,11,21,26,31,41,46,51,61,66,71,81} to 
Ati G [8, 10 19 ] and the initial state to Ati = 15. The bounds for the remaining 
entries are as before, and the remaining measurement time points were equally 
spaced. 

A plot of the functions a% , ui and 71 , 72 in the initial state and for the solution 
of problem Pq X are shown in Figure [7] A plot for the same functions with the 
same solution design as for problem Pqi after the next robustification step is 
shown in Figure |U The final design is also shown in Figure [5] A plot of the 
robustification gap Aug and the objective value of problem Pqn for each iter- 
ation N of Algorithm [T] are shown in Figure |9] A selection of design variables 
as solutions of problem Pqn is shown in Figure HHright). 

4.2 Discriminating design for two models describing signal 
sensing in dictyostelium discoideum 

The second test case is the discrimination of two models describing the chemo- 
tactic response in the amoeba Dictyostelium discoideum as presented in |33| 
using the framework presented in Section [2] and [3] The two models describe the 
adaption mechanism observed when amoebae encounter the chemoattractant 
cAMP [31], see figure [TO] For both models, a chemotaxis response regulator R 
gets activated (R*) by an activator enzyme A, when a cAMP ligand S appears. 
But the deactivating mechanism determined by the interaction with an inhibitor 
molecule / differs for both models. Both models comprise mass action kinetics 
in form of ODE. 
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Figure 5: In the left figure the robustihcation gap Ana is plotted versus the 
number of iterations TV of Algorithm [T] and in the right figure the objective value 
of problem Pgw is shown for the glycolytic design setup without the possibility 
to disturb the system. 
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Figure 6: A selection of design variables as solutions of problem Pq« for the 
glycolytic design setup without the possibility to disturb the system (left) and 
with the possibility to disturb the system (rigth) are shown. 



In model 1 the activator enzyme as well as the inhibitor enzyme are regulated 
by the external signal , which is proportional to the cAMP concentration S. 
The overall model in this case is given by, 

A\ = -k^ a Ai + k a Si 

i 1 = -k- i h + k il S 1 (50) 
i?* = -(k r A-i + k- r h)Rl + k r R T A 1: 

where fc_ a , k a , k—i, k^, k r and fc_ r are the mass action rate constants and 
Rt ■= R* + R is the total amount of the response regulator. 

In model 2 the inhibitory molecule I is activated through the indirect action of 
activator A instead of direct activation by sensing ligand binding. The overall 



20 




Figure 7: The model functions cti,a 2 and 71,72 are shown before the opti- 
mization procedure (left) and after the optimization procedure of problem Pgj 
(right) for the glycolytic design setup with the possibility to disturb the system. 
One square represents one measurement time point. 

model in this case is given by, 

A 2 = -k- a A 2 + k a S 2 

i 2 = -k-ih + k i2 A 2 (51) 
R* 2 = -(k r A 2 + k- r I 2 )R* 2 + k r R T A 2 , 

where fc_ a , k a , k—i, fcj 2 , k r and fc_ r are the mass action rate constants and 
Rt := R* + R is the total amount of the response regulator. 
For modeling details we refer to [231 ■ We have extended these systems of ordi- 
nary differential equations by an additional state corresponding to the cAMP 
ligand S with 5 = 0. By allowing species concentration perturbations c only 
to the state S we can mimic a piecewise constant control of the system by the 
cAMP ligand S. 

The experimental design parameters are the initial species concentrations of 
the four states namely, Aj, Ij, Rj, Si, the measurement time points t and the 
species concentration perturbation c with respect to S. We discard the condi- 
tion that either a measurement or a perturbation can be performed since in 
that setting by use of the perturbations c we mimic a piecewise constant input 
control S and therefore that restriction seems unnatural. Again for simplicity 
we consider the homoscedastic case with equal variances i.e. V\ = v 2 = c 2 , 
where I (Pi : P 2 ,Oi) reduces now to 

n 

l(Pi : P 2 ,Oi) = H\U) ((4 - A 2 f + (I[ - P 2 f + (i?^ - R* 2 y) . (52) 
»=i 

The parameters k_ ai k a , k-i, k^, k r , k- r and Rt are regarded as known and 
fixed, their values are given in Table [3] Parameter ki 2 is regarded as un- 
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Figure 8: The model functions ai, 0,1 and 71, 72 are shown for the same solution 
design as for problem after the next robustification step (left) and for the 
final design (right) for the glycolytic design setup with the possibility to disturb 
the system. One square represents one measurement time point. 





k a 


k-i 


k% x 




k — p 


Rt 


2.0 


3.0 


0.1 


1.0 


1.0 


1.0 


23/30 



Table 3: Parameter values for the fix values within model 1 and model 2. 

known and subject to robustification. The range of the parameter ki 2 is set 
to k l2 G [0,2]. 

The optimal design is calculated within a fixed time window with T on d = 100. 
100 equally spaced possible measurement points are defined in the initial state of 
the optimization procedure. The distance vector At between time points is sub- 
ject to design and each entry is restricted to Ai,; e [10 -7 , 10 19 ], i <G {1, 100}. 
The free perturbation vectors q, i G {11, 21, 31, 41, 51, 61, 71, 81, 91} are not re- 
stricted. The initial values are set to Cj = 0, i G {11, 21}, C31 = 0.3, c,; = —0.48, 
i G {41, 61, 81} and c, = 0.48, i G {51, 71, 91}. 

The initial species concentrations which are also subject to the experimental 
design are restricted to 57 G [0.01,0.5], Aj G [10" 7 ,1], h G [10~ 7 , 1] and 
Ri G [10~ 7 ,/]. The initial values are set to Sj = 0.2, Ai = 1.0, I = 10~ 4 and 
R = 10~ 4 . The multiple shooting intermediate variables for the species S are 
restricted to Sj G [0.01,0.5] to restrict the piecewise constant control to this in- 
terval. The parameters of the switching functions H'(ti) are chosen as a± = 5.0 
and bi = 2.5. The algorithmic settings are summarized in Table HI 
With these design conditions we start the optimization procedure twice. First 
by use of the homotopy strategy for successive problems Pq«+i with 10 homo- 
topy steps. 

Since the "discriminating power" of the experimental setup is very low in this 
case, i.e. the deviation between the two models is small, we plot the distance 
functions (Si — S2), {A\ — A2), (ii — 12) and (i?i — R2) for the initial state and 
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Figure 9: In the left figure the robustiheation gap Ana is plotted versus the 
number of iterations N of Algorithm [T] and in the right figure the objective value 
of problem Pgw is shown for the glycolytic design setup with the possibility to 
disturb the system. 




(a) model 1 (b) model 2 



Figure 10: Two models of the signal system of the Dictyostelium amoeba. 



for the solution of problem Pg^ in Figure [TTJ A plot for the same functions 
with the same solution design as for problem Pg x after the next robustification 
step is shown in Figure [121 The final design is also shown in Figure [12] 
A plot of the robustification gap A^c for each iteration TV of Algorithm [TJ is 
shown in Figure [T2] (left). A plot of the objective value of problem Pg s for each 
iteration N of Algorithm Q] is shown in Figure [H] (left). A selection of design 
variables as solutions of problem Pgw is shown in Figure IT51 (left 1 ) . 

Second we calculate the design without the homotopy strategy. We expe- 
rience huge jumps in the final objective value of problem Pgw for subsequent 
iterations ./V of Algorithm [TJ This is due to the fact that the final design of the 
former problem P§w is an infeasible starting point for the successive problem 
Pqn+i in the interior point solution strategy. First the optimizer tries to force 
the iterates back into the feasible region and afterwards the new central path 
leads to a different design. 

For this case a plot of the robustification gap Aug for each iteration N of Algo- 
rithm [TJ is shown in Figure [T3] (right). A plot of the objective value of problem 
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Optimization settings 


Integrator settings 


P 

5 


5 

kt 8 


IPOPT-tol: Step l./Step 2. 

io- 10 /io- n 


rclTol/absTol 

io- 14 /io- 14 


relTolScns/absTolSens 

io- 14 /io- 14 



Table 4: On the left hand side the optimization settings are listed comprising 
the IPOPT stopping tolerances for Step 1. and Step 2. of Algorithm [T] and on 
the right hand side the integration tolerances for the nominal trajectory and the 
first order sensitivities. We use the IPOPT option "honor.originaLbounds^io" 
for Step 1. and Step 2. of Algorithm [T] 
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Figure 11: The model variable distance functions (Si — S2), (Ai —A2), (Ii~ I2) 
and (Ri — R2) are shown before the optimization procedure (left) and after the 
optimization procedure of problem Pq X (right) for two models describing signal 
sensing in dictyostclium discoideum. One square represents one measurement 
time point. 



Pqn for each iteration N of Algorithm [T] is shown in Figure H~41 (right). A se- 
lection of design variables as solutions of problem Pgw is shown in Figure [T5] 
(right). As one can clearly see, the homotopy strategy helps to considerably 
stabilize Algorithm [TJ 



5 Conclusion 

We present a framework for the robust computation of optimal experimental de- 
signs for the purpose of model discrimination. The theoretical framework as well 
as the numerical realization by utilization of an outer approximation algorithm 
are worked out. A strategy for the numerical stabilization of the algorithm 
by use of a homotopy approach is suggested. The optimization procedure is 
successfully exemplified on two biological model systems. In our examples we 
clearly found that the homotopy approach is significantly superior to a cold 
start of successive design problems Phn+i . For the first test case, the discrimi- 



24 



0.8 1 0.8 r 1 




-0.4 ' ' ' ' ' -0.2 ' ' ' ' ' -0.4 ' ' ' ' ' -0.2 ' ' ' ' ' 

20 40 60 80 100 20 40 60 80 100 20 40 60 80 100 20 40 60 80 100 

time time time time 



Figure 12: The model variable distance functions (Si — S2), (Ai — A2), {I\ — I2) 
and — R2) are shown for the same solution design as for problem Pgj after 
the next robustification step and for the final design (right) for two models 
describing signal sensing in dictyostelium discoidcum. One square represents 
one measurement time point. 

nation of two models describing glycolytic oscillations, the outer approximation 
scheme completely fails to reach the desired accuracy 8 without homotopy strat- 
egy. For the second test case, the discrimination of two models describing signal 
sensing in dictyostelium discoidcum, the outer approximation scheme also fails 
without warmstart, however the homotopy strategy also works with only two 
homotopy steps (not presented in this paper) . This indicates the need of a step 
size strategy for reasons of efficiency which will be a next step in our work. 
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